/*
 * Copyright 2009-2019 The VOTCA Development Team (http://www.votca.org)
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */
#define BOOST_TEST_MAIN

#define BOOST_TEST_MODULE rpa_test
#include <boost/test/unit_test.hpp>
#include <votca/xtp/aobasis.h>
#include <votca/xtp/aomatrix.h>
#include <votca/xtp/logger.h>
#include <votca/xtp/orbitals.h>
#include <votca/xtp/rpa.h>
#include <votca/xtp/threecenter.h>

using namespace votca::xtp;
using namespace votca;
using namespace std;

BOOST_AUTO_TEST_SUITE(rpa_test)

BOOST_AUTO_TEST_CASE(rpa_calcenergies) {
  Logger log;
  TCMatrix_gwbse Mmn{log};
  Eigen::VectorXd eigenvals;
  RPA rpa(log, Mmn);
  rpa.configure(4, 0, 9);
  Eigen::VectorXd dftenergies = Eigen::VectorXd::Zero(10);
  dftenergies << -0.5, -0.4, -0.3, -0.2, -0.2, -0.1, 0, 0.1, 0.2, 0.3;
  Eigen::VectorXd gwenergies = Eigen::VectorXd::Zero(7);
  gwenergies << -0.15, -0.05, 0.05, 0.15, 0.45, 0.55, 0.65;
  votca::Index qpmin = 1;
  rpa.UpdateRPAInputEnergies(dftenergies, gwenergies, qpmin);
  Eigen::VectorXd rpaenergies = rpa.getRPAInputEnergies();
  Eigen::VectorXd rpaenergies_ref = Eigen::VectorXd::Zero(10);
  rpaenergies_ref << -0.85, -0.15, -0.05, 0.05, 0.15, 0.45, 0.55, 0.65, 0.75,
      0.85;
  bool e_check = rpaenergies_ref.isApprox(rpaenergies, 0.0001);

  if (!e_check) {
    cout << "energy" << endl;
    cout << rpaenergies << endl;
    cout << "energy_ref" << endl;
    cout << rpaenergies_ref << endl;
  }
  BOOST_CHECK_EQUAL(e_check, true);
}

BOOST_AUTO_TEST_CASE(rpa_full) {

  ofstream xyzfile("molecule.xyz");
  xyzfile << " 5" << endl;
  xyzfile << " methane" << endl;
  xyzfile << " C            .000000     .000000     .000000" << endl;
  xyzfile << " H            .629118     .629118     .629118" << endl;
  xyzfile << " H           -.629118    -.629118     .629118" << endl;
  xyzfile << " H            .629118    -.629118    -.629118" << endl;
  xyzfile << " H           -.629118     .629118    -.629118" << endl;
  xyzfile.close();

  ofstream basisfile("3-21G.xml");
  basisfile << "<basis name=\"3-21G\">" << endl;
  basisfile << "  <element name=\"H\">" << endl;
  basisfile << "    <shell scale=\"1.0\" type=\"S\">" << endl;
  basisfile << "      <constant decay=\"5.447178e+00\">" << endl;
  basisfile << "        <contractions factor=\"1.562850e-01\" type=\"S\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "      <constant decay=\"8.245470e-01\">" << endl;
  basisfile << "        <contractions factor=\"9.046910e-01\" type=\"S\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "    </shell>" << endl;
  basisfile << "    <shell scale=\"1.0\" type=\"S\">" << endl;
  basisfile << "      <constant decay=\"1.831920e-01\">" << endl;
  basisfile << "        <contractions factor=\"1.000000e+00\" type=\"S\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "    </shell>" << endl;
  basisfile << "  </element>" << endl;
  basisfile << "  <element name=\"C\">" << endl;
  basisfile << "    <shell scale=\"1.0\" type=\"S\">" << endl;
  basisfile << "      <constant decay=\"1.722560e+02\">" << endl;
  basisfile << "        <contractions factor=\"6.176690e-02\" type=\"S\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "      <constant decay=\"2.591090e+01\">" << endl;
  basisfile << "        <contractions factor=\"3.587940e-01\" type=\"S\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "      <constant decay=\"5.533350e+00\">" << endl;
  basisfile << "        <contractions factor=\"7.007130e-01\" type=\"S\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "    </shell>" << endl;
  basisfile << "    <shell scale=\"1.0\" type=\"SP\">" << endl;
  basisfile << "      <constant decay=\"3.664980e+00\">" << endl;
  basisfile << "        <contractions factor=\"-3.958970e-01\" type=\"S\"/>"
            << endl;
  basisfile << "        <contractions factor=\"2.364600e-01\" type=\"P\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "      <constant decay=\"7.705450e-01\">" << endl;
  basisfile << "        <contractions factor=\"1.215840e+00\" type=\"S\"/>"
            << endl;
  basisfile << "        <contractions factor=\"8.606190e-01\" type=\"P\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "    </shell>" << endl;
  basisfile << "    <shell scale=\"1.0\" type=\"SP\">" << endl;
  basisfile << "      <constant decay=\"1.958570e-01\">" << endl;
  basisfile << "        <contractions factor=\"1.000000e+00\" type=\"S\"/>"
            << endl;
  basisfile << "        <contractions factor=\"1.000000e+00\" type=\"P\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "    </shell>" << endl;
  basisfile << "  </element>" << endl;
  basisfile << "</basis>" << endl;
  basisfile.close();

  Orbitals orbitals;
  orbitals.QMAtoms().LoadFromFile("molecule.xyz");
  BasisSet basis;
  basis.Load("3-21G.xml");

  AOBasis aobasis;
  aobasis.Fill(basis, orbitals.QMAtoms());

  Eigen::VectorXd eigenvals = Eigen::VectorXd::Zero(17);
  eigenvals << 0.0468207, 0.0907801, 0.0907801, 0.104563, 0.592491, 0.663355,
      0.663355, 0.768373, 1.69292, 1.97724, 1.97724, 2.50877, 2.98732, 3.4418,
      3.4418, 4.81084, 17.1838;

  Eigen::MatrixXd eigenvectors = Eigen::MatrixXd::Zero(17, 17);
  eigenvectors << 0.0185815, 2.9133e-17, 8.49354e-17, -0.00312916, 0.0420075,
      1.11356e-16, 1.85886e-17, -0.0334732, 0.0485113, -8.71556e-18,
      -3.79994e-17, -0.0346485, -0.0248392, -3.32286e-22, -4.62643e-17,
      0.0144472, 0.996183, 0.166534, 2.80578e-16, 7.85515e-16, -0.0299557,
      0.409156, 1.10045e-15, 1.67666e-16, -0.336895, 0.608646, -1.00343e-16,
      -5.04622e-16, -0.437568, -0.299751, -2.77682e-17, -6.06812e-16, 0.176594,
      -0.0866677, 0.0010572, -0.0210402, -0.0345975, 0.035778, 0.0611836,
      0.0374747, -0.154443, 0.0892921, 0.0842611, -0.35309, -0.0759572,
      0.278374, -0.409082, -0.64367, 0.308248, -0.261525, -0.000315534,
      -0.0010572, -0.0404824, 0.000922593, -0.035778, -0.0611836, -0.115015,
      -0.109676, -0.0892921, -0.0842611, -0.242326, 0.267806, -0.278374,
      0.409082, -0.0548848, 0.711558, 0.261525, 0.000315534, 0.0010572,
      -0.0194422, 0.0355201, 0.035778, 0.0611836, -0.152489, 0.0447677,
      0.0892921, 0.0842611, 0.110764, 0.343764, 0.278374, -0.409082, 0.588785,
      0.403311, -0.261525, -0.000315534, -0.823783, -9.8891e-16, -3.34692e-15,
      0.103497, -0.277613, -5.51463e-16, 4.95594e-17, 0.163544, 0.121215,
      -7.23985e-17, -2.63149e-16, -0.259891, -0.284396, -1.74149e-16,
      -6.65818e-16, 0.208987, 0.00782842, -0.0333718, 0.22696, 0.373203,
      -0.337332, -0.251625, -0.144131, 0.594004, -0.329076, -0.0456626, 0.18588,
      0.0399869, -0.0631275, -0.0704844, -0.231899, 0.111054, -0.189161,
      0.000129868, 0.0333718, 0.436683, -0.009952, 0.337332, 0.251625, 0.442357,
      0.421824, 0.329076, 0.0456626, 0.12757, -0.140984, 0.0631275, 0.0704844,
      -0.0197737, 0.256357, 0.189161, -0.000129868, -0.0333718, 0.209723,
      -0.383155, -0.337332, -0.251625, 0.586489, -0.172181, -0.329076,
      -0.0456626, -0.0583106, -0.180971, -0.0631275, -0.0704844, 0.212125,
      0.145303, -0.189161, 0.000129868, -0.00177478, 0.0553645, -0.00126176,
      -0.0164247, 0.23154, -0.262519, -0.250334, -0.0135392, -0.429472, 0.45567,
      -0.503583, -0.223493, -0.211802, -0.020461, 0.265268, 0.0023362,
      -0.00241145, 0.294363, -0.686239, 0.0156394, 0.204055, -0.360136,
      0.267096, 0.254698, 0.074687, -0.0228668, 0.132236, -0.14614, -0.174986,
      -0.185046, -0.0109958, 0.142556, 0.0661743, 0.0022999, -0.00177478,
      -0.0265895, 0.0485779, -0.0164247, 0.23154, 0.348055, -0.102182,
      -0.0135392, -0.429472, 0.208281, 0.646413, -0.223493, -0.211802,
      -0.219498, -0.150354, 0.0023362, -0.00241145, 0.294363, 0.329576,
      -0.60212, 0.204055, -0.360136, -0.354123, 0.103963, 0.074687, -0.0228668,
      0.0604434, 0.18759, -0.174986, -0.185046, -0.117959, -0.0808008,
      0.0661743, 0.0022999, -0.00177478, -0.028775, -0.0473162, -0.0164247,
      0.23154, -0.0855356, 0.352515, -0.0135392, -0.429472, -0.663951, -0.14283,
      -0.223493, -0.211802, 0.239959, -0.114914, 0.0023362, -0.00241145,
      0.294363, 0.356664, 0.586481, 0.204055, -0.360136, 0.0870267, -0.358661,
      0.074687, -0.0228668, -0.192679, -0.0414494, -0.174986, -0.185046,
      0.128955, -0.0617554, 0.0661743, 0.0022999, 0.00741062, -3.87173e-16,
      -4.31863e-16, 0.0468488, -0.0476991, 7.27357e-16, 1.23654e-15, -0.43422,
      -0.159247, -4.34945e-17, 1.2743e-16, 0.503528, -0.228856, -7.97629e-17,
      -2.53026e-16, 0.689669, -0.00301027, 0.173046, 7.91486e-15, 8.39419e-15,
      -0.717804, -0.0195249, -1.10754e-15, -1.66789e-15, 0.551371, 0.0684292,
      4.15572e-17, -1.84233e-16, 0.0105378, -0.148396, -1.63792e-16,
      -4.6499e-16, 0.351571, 0.00210309;

  Logger log;
  TCMatrix_gwbse Mmn{log};
  Mmn.Initialize(aobasis.AOBasisSize(), 0, 16, 0, 16);
  Mmn.Fill(aobasis, aobasis, eigenvectors);

  RPA rpa(log, Mmn);
  rpa.configure(4, 0, 16);
  rpa.setRPAInputEnergies(eigenvals);
  Eigen::MatrixXd e_i = rpa.calculate_epsilon_i(0.5);

  Eigen::MatrixXd i_ref = Eigen::MatrixXd::Zero(17, 17);
  i_ref << 1.00063, -0.000738023, -0.000847951, 0.000847951, -0.000847951,
      -0.00398178, -0.00230038, 0.00230038, -0.00230038, -0.000509516,
      -0.00507621, -0.000509516, -0.00507621, -0.000509516, -0.00507621,
      -0.000383038, -0.00312618, -0.000738023, 1.00263, 0.00161324, -0.00161324,
      0.00161324, 0.00911444, 0.00415952, -0.00415952, 0.00415952, 0.000514703,
      0.00743779, 0.000514703, 0.00743779, 0.000514703, 0.00743779, 0.000660954,
      0.00451844, -0.000847951, 0.00161324, 1.00392, -0.00190885, 0.00190885,
      0.00907072, 0.0108943, -0.00558202, 0.00558202, 0.00178259, 0.0142578,
      0.00178259, 0.0142578, 0.00103807, 0.012368, 0.00152356, 0.00935261,
      0.000847951, -0.00161324, -0.00190885, 1.00392, -0.00190885, -0.00907072,
      -0.00558202, 0.0108943, -0.00558202, -0.00103807, -0.012368, -0.00178259,
      -0.0142578, -0.00178259, -0.0142578, -0.00152356, -0.00935261,
      -0.000847951, 0.00161324, 0.00190885, -0.00190885, 1.00392, 0.00907072,
      0.00558202, -0.00558202, 0.0108943, 0.00178259, 0.0142578, 0.00103807,
      0.012368, 0.00178259, 0.0142578, 0.00152356, 0.00935261, -0.00398178,
      0.00911444, 0.00907072, -0.00907072, 0.00907072, 1.06083, 0.0317799,
      -0.0317799, 0.0317799, 0.00979314, 0.0779294, 0.00979314, 0.0779294,
      0.00979314, 0.0779294, 0.00426666, 0.0435389, -0.00230038, 0.00415952,
      0.0108943, -0.00558202, 0.00558202, 0.0317799, 1.06312, -0.0152029,
      0.0152029, 0.0130855, 0.0740544, 0.0130855, 0.0740544, -0.00319113,
      0.0230598, 0.0056834, 0.0299224, 0.00230038, -0.00415952, -0.00558202,
      0.0108943, -0.00558202, -0.0317799, -0.0152029, 1.06312, -0.0152029,
      0.00319113, -0.0230598, -0.0130855, -0.0740544, -0.0130855, -0.0740544,
      -0.0056834, -0.0299224, -0.00230038, 0.00415952, 0.00558202, -0.00558202,
      0.0108943, 0.0317799, 0.0152029, -0.0152029, 1.06312, 0.0130855,
      0.0740544, -0.00319113, 0.0230598, 0.0130855, 0.0740544, 0.0056834,
      0.0299224, -0.000509516, 0.000514703, 0.00178259, -0.00103807, 0.00178259,
      0.00979314, 0.0130855, 0.00319113, 0.0130855, 1.01191, 0.0403377,
      0.000522925, 0.0147762, 0.000522925, 0.0147762, 0.0012644, 0.0119637,
      -0.00507621, 0.00743779, 0.0142578, -0.012368, 0.0142578, 0.0779294,
      0.0740544, -0.0230598, 0.0740544, 0.0403377, 1.19801, 0.0147762, 0.12575,
      0.0147762, 0.12575, 0.00894661, 0.0796677, -0.000509516, 0.000514703,
      0.00178259, -0.00178259, 0.00103807, 0.00979314, 0.0130855, -0.0130855,
      -0.00319113, 0.000522925, 0.0147762, 1.01191, 0.0403377, 0.000522925,
      0.0147762, 0.0012644, 0.0119637, -0.00507621, 0.00743779, 0.0142578,
      -0.0142578, 0.012368, 0.0779294, 0.0740544, -0.0740544, 0.0230598,
      0.0147762, 0.12575, 0.0403377, 1.19801, 0.0147762, 0.12575, 0.00894661,
      0.0796677, -0.000509516, 0.000514703, 0.00103807, -0.00178259, 0.00178259,
      0.00979314, -0.00319113, -0.0130855, 0.0130855, 0.000522925, 0.0147762,
      0.000522925, 0.0147762, 1.01191, 0.0403377, 0.0012644, 0.0119637,
      -0.00507621, 0.00743779, 0.012368, -0.0142578, 0.0142578, 0.0779294,
      0.0230598, -0.0740544, 0.0740544, 0.0147762, 0.12575, 0.0147762, 0.12575,
      0.0403377, 1.19801, 0.00894661, 0.0796677, -0.000383038, 0.000660954,
      0.00152356, -0.00152356, 0.00152356, 0.00426666, 0.0056834, -0.0056834,
      0.0056834, 0.0012644, 0.00894661, 0.0012644, 0.00894661, 0.0012644,
      0.00894661, 1.00208, 0.00789381, -0.00312618, 0.00451844, 0.00935261,
      -0.00935261, 0.00935261, 0.0435389, 0.0299224, -0.0299224, 0.0299224,
      0.0119637, 0.0796677, 0.0119637, 0.0796677, 0.0119637, 0.0796677,
      0.00789381, 1.06135;

  bool i_check = i_ref.isApprox(e_i, 0.0001);

  if (!i_check) {
    cout << "Epsilon_i" << endl;
    cout << e_i << endl;
    cout << "Epsilon_i_ref" << endl;
    cout << i_ref << endl;
  }
  BOOST_CHECK_EQUAL(i_check, 1);

  Eigen::MatrixXd e_r = rpa.calculate_epsilon_r(0.0);

  Eigen::MatrixXd r_ref = Eigen::MatrixXd::Zero(17, 17);
  r_ref << 1.0006587695, -0.00078799598102, -0.00088922577334, 0.00088922577334,
      -0.00088922577334, -0.0042169065833, -0.002336091509, 0.002336091509,
      -0.002336091509, -0.00050662218169, -0.0052057890458, -0.00050662218169,
      -0.0052057890458, -0.00050662218169, -0.0052057890458, -0.00037564339451,
      -0.0032058806556, -0.00078799598102, 1.0028668699, 0.0017910236524,
      -0.0017910236524, 0.0017910236524, 0.01014801061, 0.0051408879134,
      -0.0051408879134, 0.0051408879134, 0.00036877585664, 0.0073965347948,
      0.00036877585664, 0.0073965347948, 0.00036877585664, 0.0073965347948,
      0.00094029956345, 0.0045022793469, -0.00088922577334, 0.0017910236524,
      1.0067257854, -0.0015465163998, 0.0015465163998, 0.0094511460318,
      0.045585007947, 0.0041655130822, -0.0041655130822, 0.0034824005183,
      0.028928879976, 0.0034824005183, 0.028928879976, -0.003242803598,
      -0.014354267669, 0.0031137836848, 0.0099736516652, 0.00088922577334,
      -0.0017910236524, -0.0015465163998, 1.0067257854, -0.0015465163998,
      -0.0094511460318, 0.0041655130822, 0.045585007947, 0.0041655130822,
      0.003242803598, 0.014354267669, -0.0034824005183, -0.028928879976,
      -0.0034824005183, -0.028928879976, -0.0031137836848, -0.0099736516652,
      -0.00088922577334, 0.0017910236524, 0.0015465163998, -0.0015465163998,
      1.0067257854, 0.0094511460318, -0.0041655130822, 0.0041655130822,
      0.045585007947, 0.0034824005183, 0.028928879976, -0.003242803598,
      -0.014354267669, 0.0034824005183, 0.028928879976, 0.0031137836848,
      0.0099736516652, -0.0042169065833, 0.01014801061, 0.0094511460318,
      -0.0094511460318, 0.0094511460318, 1.0683943743, 0.033253778547,
      -0.033253778547, 0.033253778547, 0.0098928173624, 0.080323544914,
      0.0098928173624, 0.080323544914, 0.0098928173624, 0.080323544914,
      0.0042490675042, 0.044368325677, -0.002336091509, 0.0051408879134,
      0.045585007947, 0.0041655130822, -0.0041655130822, 0.033253778547,
      1.5522819346, 0.1630437813, -0.1630437813, 0.040951984157, 0.29148484041,
      0.040951984157, 0.29148484041, -0.065920061912, -0.39062789159,
      0.019262830405, 0.030434350374, 0.002336091509, -0.0051408879134,
      0.0041655130822, 0.045585007947, 0.0041655130822, -0.033253778547,
      0.1630437813, 1.5522819346, 0.1630437813, 0.065920061912, 0.39062789159,
      -0.040951984157, -0.29148484041, -0.040951984157, -0.29148484041,
      -0.019262830405, -0.030434350374, -0.002336091509, 0.0051408879134,
      -0.0041655130822, 0.0041655130822, 0.045585007947, 0.033253778547,
      -0.1630437813, 0.1630437813, 1.5522819346, 0.040951984157, 0.29148484041,
      -0.065920061912, -0.39062789159, 0.040951984157, 0.29148484041,
      0.019262830405, 0.030434350374, -0.00050662218169, 0.00036877585664,
      0.0034824005183, 0.003242803598, 0.0034824005183, 0.0098928173624,
      0.040951984157, 0.065920061912, 0.040951984157, 1.0212195523,
      0.09982358437, -0.0032454628838, -0.012437519198, -0.0032454628838,
      -0.012437519198, 0.00041150551217, 0.012585676779, -0.0052057890458,
      0.0073965347948, 0.028928879976, 0.014354267669, 0.028928879976,
      0.080323544914, 0.29148484041, 0.39062789159, 0.29148484041,
      0.09982358437, 1.6097216877, -0.012437519198, -0.062989343739,
      -0.012437519198, -0.062989343739, 0.010200010299, 0.082551786343,
      -0.00050662218169, 0.00036877585664, 0.0034824005183, -0.0034824005183,
      -0.003242803598, 0.0098928173624, 0.040951984157, -0.040951984157,
      -0.065920061912, -0.0032454628838, -0.012437519198, 1.0212195523,
      0.09982358437, -0.0032454628838, -0.012437519198, 0.00041150551217,
      0.012585676779, -0.0052057890458, 0.0073965347948, 0.028928879976,
      -0.028928879976, -0.014354267669, 0.080323544914, 0.29148484041,
      -0.29148484041, -0.39062789159, -0.012437519198, -0.062989343739,
      0.09982358437, 1.6097216877, -0.012437519198, -0.062989343739,
      0.010200010299, 0.082551786343, -0.00050662218169, 0.00036877585664,
      -0.003242803598, -0.0034824005183, 0.0034824005183, 0.0098928173624,
      -0.065920061912, -0.040951984157, 0.040951984157, -0.0032454628838,
      -0.012437519198, -0.0032454628838, -0.012437519198, 1.0212195523,
      0.09982358437, 0.00041150551217, 0.012585676779, -0.0052057890458,
      0.0073965347948, -0.014354267669, -0.028928879976, 0.028928879976,
      0.080323544914, -0.39062789159, -0.29148484041, 0.29148484041,
      -0.012437519198, -0.062989343739, -0.012437519198, -0.062989343739,
      0.09982358437, 1.6097216877, 0.010200010299, 0.082551786343,
      -0.00037564339451, 0.00094029956345, 0.0031137836848, -0.0031137836848,
      0.0031137836848, 0.0042490675042, 0.019262830405, -0.019262830405,
      0.019262830405, 0.00041150551217, 0.010200010299, 0.00041150551217,
      0.010200010299, 0.00041150551217, 0.010200010299, 1.0065147698,
      0.0087814549722, -0.0032058806556, 0.0045022793469, 0.0099736516652,
      -0.0099736516652, 0.0099736516652, 0.044368325677, 0.030434350374,
      -0.030434350374, 0.030434350374, 0.012585676779, 0.082551786343,
      0.012585676779, 0.082551786343, 0.012585676779, 0.082551786343,
      0.0087814549722, 1.0692401254;

  bool r_check = r_ref.isApprox(e_r, 0.0001);

  if (!r_check) {
    cout << "Epsilon_r" << endl;
    cout << e_r << endl;
    cout << "Epsilon_r_ref" << endl;
    cout << r_ref << endl;
  }

  BOOST_CHECK_EQUAL(r_check, 1);
}

BOOST_AUTO_TEST_SUITE_END()
